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Simulations  of  Cavity-stabilized  Flames  in  Supersonic  Flows 
Using  Reduced  Chemical  Kinetic  Mechanisms 
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Tianfeng  Lu*  and  Chung  K.  Law§ 

Department  of  Mechanical  and  Aerospace  Engineering,  Princeton  University,  Princeton,  NJ  08544 

The  VULCAN  CFD  code  integrated  with  a  reduced  chemical  kinetic  mechanism  was 
applied  to  simulate  cavity-stabilized  ethylene-air  flames  and  to  predict  flame  stability  limits 
in  supersonic  flows  based  on  an  experimental  study.  A  15-step  reduced  kinetic  mechanism 
for  ethylene  was  systematically  developed  through  skeletal  reduction  with  a  directed  relation 
graph  and  time  scale  reduction  based  on  quasi-steady  state  assumptions.  The  accuracy  of 
the  reduced  kinetic  mechanism  and  its  implementation  in  the  VULCAN  code  were 
demonstrated  in  an  auto-ignition  problem  with  a  range  of  parameters.  3D  simulations  were 
then  carried  out  for  cavity- stabilized  flames  at  different  fuel  flowrates  and  turbulent 
Schmidt  numbers.  For  comparison  with  the  performance  of  the  present  reduced 
mechanism,  a  3-  and  a  10-step  global  kinetic  model  were  applied  to  simulate  the  same  cavity 
combustor,  and  the  results  show  that  the  15-step  reduced  model  predicts  experimental 
results  much  better  than  the  3-  and  10-step  models.  The  importance  of  including  accurate 
chemical  kinetics  in  CFD  simulations  is  therefore  demonstrated. 
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immediate  error 
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global  equivalence  ratio 
reaction  rate 


I.  Introduction 

The  U.S.  Air  Force  (USAF)  has  been  the  supporting  the  development  of  hydrocarbon-fueled  scramjet  engines 
via  the  HyTech  Program  since  1995. 1  A  key  issue  that  has  been  addressed  in  this  program  is  the  fuel  and  air  mixing 
and  flameholding.  The  mixing  of  fuel  and  air  in  an  airbreathing  engine  becomes  increasingly  inefficient  and 
flameholding  becomes  increasingly  difficult  at  high  velocities,  due  to  the  reduced  flow  residence  time  in  the 
combustors  and  the  flow  compressibility,  which  adversely  affects  the  mixing  rate.  Wall  cavities  are  outstanding 
flameholders  for  scramjet  combustors  because  they  provide  sufficient  residence  time  for  the  flameholding  of 
hydrocarbon  fuels  in  a  supersonic  stream  and  also  they  involve  low  total  pressure  losses.  In  the  last  decade,  the  Air 
Force  Research  Laboratory's  Propulsion  Directorate,  Aerospace  Propulsion  Division  (AFRL/PRA)  has  been 
developing  several  modern  direct- connect  test  facilities  involving  different  types  of  cavity  flameholders  for 
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hydrocarbon- fueled  scramjet  combustors.  In  an  effort  to  develop  robust  flameholding  schemes,  flame  stability 
studies  have  been  receiving  special  attention  from  these  facilities.  For  example,  Mathur  et  al.2  demonstrated  the 
ability  of  cavity  flameholders  to  stabilize  hydrocarbon  flames  at  flight  Mach  numbers  between  4  and  6  over  a  wide 
range  of  combustor  fuel-air  equivalence  ratios.  Gruber  et  al.3  provided  limited  data  on  the  operating  range  of 
cavity- stabilized  flames.  Rasmussen,  et  al.4  obtained  stability  limits  in  directly  fueled  cavity  flameholders  for 
different  air  mass  flowrates,  cavity  geometries,  fuel  injection  schemes,  Mach  numbers,  and  fuel  types.  Most  of 
these  fundamental  researches  used  ethylene  as  a  fuel,  because  it  has  short  ignition  delays,  good  energy  yields,  and  is 
a  key  component  of  cracked  endothermic  fuels  such  as  JP7.  While  most  of  the  investigations  focused  on 
experimental  studies,  there  has  been  little  effort  in  numerical  studies  on  the  stability  of  hydrocarbon-fueled  flames 
in  cavity  flameholders  in  supersonic  airflows,  primarily  due  to  the  inherent  complexities  in  physics. 

Flame  stability  limits  are  characterized  with  either  the  fuel  lean  or  fuel  rich  blowout.  Flame  blowout  is  a  highly 
transient  phenomenon  and  is  very  sensitive  to  aerodynamic  and  thermodynamic  conditions.  For  instance,  blowout 
is  generally  abrupt  as  a  result  of  the  Arrhenius  response  of  chemical  kinetics  in  competition  with  the  more  gradual 
heat  loss  through  fluid  mixing.  Therefore,  it  is  crucial  to  incorporate  accurate  finite-rate  chemistry  in  the  form  of 
either  detailed  or  reduced  mechanisms  for  accurate  prediction  of  blowout,  in  that  even  a  small  or  moderate 
uncertainty  in  chemical  kinetic  rates  could  lead  to  gross  errors  in  the  simulation  results.  However,  detailed  reaction 
mechanisms  are  typically  too  computationally  demanding  due  to  the  large  number  of  species  and  reactions  required, 
as  well  as  the  stiffness  induced  by  the  vastly  different  time  scales,  which  may  span  many  orders  of  magnitudes, 
between  the  major  species  and  the  highly  reactive  radicals.  As  such,  detailed  chemistry  is  rarely  included  in 
complex  3D  turbulent  flame  simulations  even  for  the  simplest  hydrocarbons,  such  as  ethylene,  and  even  by  using 
the  latest  supercomputers.  To  overcome  this  difficulty,  extensive  efforts  have  been  made  to  develop  reduced 
reaction  mechanisms  of  smaller  sizes  and  reduced  stiffness,  while  still  being  sufficiently  accurate.  Reduced 
mechanisms  are  typically  derived  from  validated  detailed  mechanisms  at  two  levels,  namely  skeletal  reduction, 
which  eliminates  unimportant  reactions  and  species,  and  time  scale  reduction,  which  involves  partial  equilibrium 
and  quasi-steady  state  (QSS)  assumptions. 

The  objective  of  the  current  study  is  to  simulate  cavity- stabilized  flames  and  to  predict  flame  stability  limits  in 
supersonic  flows  using  reduced  chemical  kinetic  mechanisms  based  on  an  experimental  study  with  the  AFRL/PRA 
direct-connect  cavity  combustor  facility.  The  experimental  study  focused  on  the  stability  of  ethylene-fueled  flames 
and  the  results  were  reported  in  Ref.  4.  The  reduced  mechanism  in  the  present  study  was  developed  using  the 
directed  relation  graph  (DRG)  method  for  skeletal  reduction  and  the  computational  singular  perturbation  (CSP) 
theory  for  QSS  based  time  scale  reduction.5  This  15-step  mechanism  was  integrated  into  VULCAN,6  which  is  a 
specialized  CFD  code  for  the  analysis  of  high  speed  turbulent  reacting  flows  and  has  been  widely  adopted  for 
simulations  of  scramjet  engines.  A  brief  review  of  the  reduction  methods  and  the  resulting  mechanism  will  be 
provided  in  the  following  section,  followed  by  the  description  of  VULCAN.  The  VULCAN  code  integrated  with 
the  reduced  mechanism  is  first  validated  in  an  auto-ignition  problem,  and  is  then  applied  to  simulate  the  AFRL/PRA 
cavity  combustor  at  a  selected  test  condition  with  different  fuel  flowrates  and  turbulent  Schmidt  numbers.  In  the 
end,  to  demonstrate  the  importance  of  using  an  accurate  reduced  kinetic  model  in  the  current  simulation,  two  global 
kinetic  models,  with  3  and  10  steps  respectively,  will  be  also  applied  to  the  same  cavity  combustor,  and  the 
discrepancies  among  the  results  from  the  different  models  will  be  discussed. 


II.  Methodology  for  Mechanism  Reduction 


A.  Skeletal  reduction  using  directed  relation  graph 


The  skeletal  reduction  is  to  identify  and  eliminate  unimportant  reactions  or  species.  The  elimination  of  species  is 
of  primary  interest  for  most  combustion  simulations  involving  a  Jacobian  matrix  which  scales  quadratically  with  the 
number  of  species.  The  method  of  DRG  was  developed  specifically  for  direct  elimination  of  unimportant  species 
by  resolving  species  couplings  through  graph  searching.  In  DRG,  direct  species  coupling  was  defined  by  the  relative 
immediate  error  induced  to  one  species,  say  A,  due  to  the  elimination  of  another,  say  B: 
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(1) 


where  the  subscripts  A  and  B  indicate  the  species  identity,  i  is  the  ith  reaction,  Va,i  the  stoichiometric  coefficient  of 
species  A  in  the  ith  reaction,  and  CO  the  reaction  rate.  It  is  seen  that  the  denominator  consists  of  the  contributions  to 
species  A  from  each  reaction,  and  the  numerator  consists  of  all  the  terms  involving  species  B  in  the  denominator. 
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Therefore  vAb  roughly  describes  the  importance  of  species  B  to  A.  If  the  relative  error  rAB  is  not  small  compared  to  a 
threshold  value  8,  the  removal  of  species  B  from  the  skeletal  mechanism  immediately  induces  a  non-negligible  error 
in  A,  and  if  A  is  to  be  retained  in  the  skeletal  mechanism,  B  should  also  be  retained  for  accurate  evaluation  of  the 
reaction  rate  of  A.  This  requirement  of  species  A  to  B  is  denoted  as  A  — »  B. 

Species  can  also  be  coupled  indirectly  through  direct  coupling  with  intermediate  species.  For  example,  if  A  — » 
B  and  B  — »  C,  then  A  requires  C  indirectly.  The  set  of  species  required  either  directly  or  indirectly  by  species  A  is 
defined  as  the  dependent  set  of  A,  and  it  is  readily  seen  that  if  species  A  is  to  be  retained  in  the  skeletal  mechanisms, 
so  should  be  all  the  species  in  its  dependent  set. 

A  sample  relation  of  DRG  is  shown  in  Fig.  1  in  which  vertices  A,  B  . . .  F 
correspond  to  the  respective  species.  The  arrow  shows  the  direction  of 
dependence  of  one  species  on  another,  with  its  width  indicating  the  strength 
of  the  dependence.  The  starting  vertex  A  is  a  major  species  in  the 
mechanism.  The  skeletal  species  in  Fig.  1  correspond  to  all  the  vertices 
reachable  from  the  starting  vertex.  Graph  searching  methods  such  as  deep 
first  search  can  thus  be  exploited  to  efficiently  find  all  the  vertices  reachable 
from  the  starting  one.  Species  A,  B  and  D  are  therefore  identified  as  the 
species  of  the  skeletal  mechanism  in  the  current  example. 

To  obtain  a  skeletal  mechanism  over  a  sufficiently  wide  range  of  Figure  1.  A  DRG  showing  typical 
parameters,  such  as  pressure,  temperature,  equivalence  ratio,  and  residence  relations  of  the  species, 

time,  a  set  of  reaction  states  are  sampled  from  typical  applications,  such  as 

perfectly  stirred  reactor  (PSR)  and  auto-ignition.  A  local  DRG  reduction  can  be  performed  for  each  of  these 
reaction  states,  and  the  union  of  the  species  sets  of  all  the  local  skeletal  mechanisms  yields  the  species  set  of  the 
global  skeletal  mechanism.  Strictly  speaking,  the  skeletal  mechanism  developed  using  DRG  should  only  be  applied 
to  those  applications  and  parametric  spaces  based  on  which  it  is  developed.  However,  it  is  often  the  case  that  a 
skeletal  mechanism  developed  from  simple  applications  can  be  applied  to  more  complicated  simulations  and  still 
has  satisfactory  accuracy.7 

B.  Further  reduction  based  on  QSS  assumptions 

Since  a  skeletal  mechanism  contains  only  elementary  reactions,  it  is  applicable  in  the  same  manner  as  the 
detailed  mechanism.  However,  the  size  of  a  skeletal  mechanism  is  usually  still  too  large  for  many  demanding 
simulations.  A  further  reduction  can  be  accomplished  by  eliminating  species  with  short  time  scales  by 
approximating  them  to  be  in  a  steady  state.  The  removal  of  the  short  time  scales  effectively  reduces  the  stiffness  of 
the  system  as  well  as  the  number  of  differential  equations.  This  combined  benefit  usually  dominates  over  the  costs 
incurred  by  the  internal  algebraic  iteration  in  solving  the  concentrations  of  the  QSS  species,  and  consequently 
results  in  an  overall  time  saving  in  using  QSS-based  reduced  mechanisms. 

The  short  time  scales  can  be  identified  by  several  methods  and  the  CSP  method  was  used  in  the  present  study. 
In  this  method,  the  time  scales  of  independent  modes  are  summed  to  an  indicator  weighed  by  sufficiently  large 
radical  pointers,  the  reciprocal  of  which  can  be  defined  as  the  time  scale  of  the  species.  A  critical  time  scale  is 
selected  to  normalize  the  time  scale  of  species,  such  that  the  time  scale  of  the  major  species  controlling  the  overall 
reaction  rate  of  the  system  is  normalized  to  unity.  Species  are  considered  to  be  QSS  if  their  normalized  time  scales 
are  smaller  than  a  specified  critical  value.  The  selection  of  the  critical  time  scale  depends  on  the  system.  The 
extinction  time  for  the  PSR,  and  the  auto-ignition  time,  can  be  respectively  used  as  the  characteristic  time  for 
extinction  and  ignition  problems. 

C.  Mechanism  reduction  for  ethylene 

By  applying  the  above  methodology  a  reduced  chemical  kinetic  mechanism  was  created  for  ethylene.  The 
detailed  ethylene  oxidation  mechanism  used  for  reduction  is  from  Qin  et  al.8  and  it  consists  of  70  species  and  463 
elementary  reactions.  For  skeletal  reduction,  PSR  and  auto-ignition  were  used  as  the  data  source  for  reduction  with 
pressure  ranging  from  0.1  ~  30  atm,  equivalence  ratio  from  0.7  ~  1.3,  initial  temperature  of  300  K  for  PSR  and  1000 
~  1800  K  for  auto-ignition,  and  residence  time  covering  the  entire  ignition  and  extinction  ranges.  Data  points  were 
densely  sampled  over  the  above  parametric  space  to  ensure  that  all  specific  relations  among  the  species  can  be 
represented  by  the  sample  set.  During  the  graph  searching,  ethylene  was  selected  as  the  starting  species.  By 
choosing  the  threshold  accuracy  value  8  of  0.16,  a  skeletal  mechanism  consisting  of  33  species  was  obtained.  The 
species  constituting  this  skeletal  mechanism  are:  H2,  H,  0 ,  02,  OH,  H20,  H02,  H202,  C,  CH,  CH2,  CH2*,  CH3,  CH4, 
CO,  C02,  HCO,  CH20,  CH2OH,  CHsO,  CHsOH,  C2H2,  C2H3,  C2H4,  C2H5,  HCCO,  CH2CO,  CH2CHO,  n-C3H7,  C3H6, 
a-C3H5,  Ar,  and  N2.  By  eliminating  all  elementary  reactions  in  the  detailed  mechanism  which  has  any  reactants  or 
products  not  included  in  the  above  set  of  species,  205  elementary  reactions  were  retained  for  the  skeletal 
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mechanism. 

The  33-  species  skeletal  mechanism  developed  above  was  further  reduced  by  using  the  CSP  to  identify  the  QSS 
species.  By  again  using  PSR  and  auto-ignition  as  the  data  bases,  the  time  scale  for  each  species  was  calculated  for 
each  data  point,  with  their  characteristic  times  being  those  of  the  extinction  and  ignition  time  respectively.  The 
longest  normalized  time  for  each  species  was  taken  as  the  characteristic  time  scale  for  this  species  because  it 
describes  approximately  its  worst  case  characteristic  time.  A  species  is  assumed  to  be  a  QSS  species  if  its  worst- 
case  normalized  time  scale  is  shorter  than  the  specified  threshold  value  8  of  0.025.  The  QSS  species  identified  were 
then  removed  from  the  skeletal  mechanism  and  an  internal  algebraic  loop  was  used  to  solve  the  concentrations  of 
the  QSS  species.  The  final  reduced  mechanism  consists  of  19  species  or  15  global  steps.  These  19  species  are:  H2 , 
H,  O,  02,  OH,  H20,  H02,  H202,  CH3,  CH4,  CO,  C02,  CH20,  C2H2,  C2H4,  C2H6,  CH2CO,  C3H6,  and  N2. 


III.  VULCAN  CFD  Code 

In  this  study  the  CFD  code  in  which  the  above  reduced  kinetic  mechanism  was  implemented  is  based  on  the 
VULCAN  code.  In  contrast  to  most  CFD  codes,  VULCAN  has  been  developed  from  its  inception  to  be  an  analysis 
tool  for  high-speed  turbulent  reacting  flows.  The  VULCAN  flow  solver  is  presently  being  used  in  the  USAF 
HyTech  program  as  well  as  NASA  hypersonic  vehicle  programs.  In  addition,  VULCAN  is  also  used  by  both 
industry  and  academic  institutions.  Ongoing  development  and  maintenance  of  VULCAN  is  being  performed  by 
NASA  Langley.  Its  development  is  primarily  driven  by  application  to  high-speed  propulsion  devices. 

VULCAN  is  a  density-based  CFD  code  applicable  to  complicated  2D  and  3D  geometries  by  using  multiblock 
structured  grid.  Grid  lines  across  a  block  interface  can  be  continuous  or  discontinuous.  The  code  solves  the 
Reynolds-averaged  Navier  Stokes  (RANS)  equations  appropriate  for  calorically  or  thermally  perfect  gases  with  a 
cell-centered  finite  volume  scheme.  The  equation  set  can  be  integrated  in  a  fully  elliptic  or  space-marched  manner. 
The  in  viscid  fluxes  can  be  evaluated  with  central  differences,  Roe’s  flux  difference  method,  or  a  low-diffusion  flux 
vector  split  scheme.  Several  flux  limiters  are  provided  to  ensure  total  variation  diminishing.  A  variety  of  two- 
equation  turbulence  models  are  available,  along  with  one-equation  and  explicit  algebraic  Reynolds  stress  models. 
Assumed  PDF  options  exist  for  modeling  turbulence-chemistry  interactions.  Chemically  reacting  flows  can  be 
modeled  with  general  finite-rate  kinetics  models.  The  code  also  contains  full  multi-grid  capabilities,  allowing  rapid 
convergence  for  steady-state  problems.  For  parallel  computation,  the  MPI  message  passing  library  is  used  as  a 
communication  software  in  the  code. 

In  the  original  VULCAN  code  the  chemical  kinetic  models  must  be  based  on  the  standard  Arrhenius  equation 
and  the  reaction  rate  depends  only  on  temperature.  In  reality,  reaction  rates  may  depend  on  pressure  as  well  as 
temperature.  In  the  present  reduced  kinetic  mechanism  of  ethylene,  either  the  Lindemann  or  Troe  expression  is 
applied  to  account  for  the  pressure  dependence  in  some  reactions.  In  order  to  accommodate  the  pressure  dependent 
reaction  rate,  the  VULCAN  code  was  modified  by  introducing  a  subroutine  as  a  driver.  This  driver  is  called  to 
provide  flow  information  such  as  temperature,  pressure,  species  concentrations,  etc.  every  time  before  a  reaction 
rate  is  computed.  Introducing  this  driver  leads  to  little  change  in  the  source  code. 

Turbulence  modeling  is  one  of  the  largest  challenges  in  numerical  analysis  of  turbulent  reacting  flows.  In  this 
study  the  turbulence  model  was  based  on  the  Menter  turbulence  model.  The  Menter  model  is  essentially  the 
standard  high-Reynolds-number  form  of  the  Wilcox  k-co  model  near  solid  surfaces,  but  it  smoothly  transitions  to  a 
standard  Jones-Launder  k-8  model  near  the  outer  portion  of  the  boundary  layer  and  in  regions  of  free  shear.  To 
relax  the  grid  requirements  near  solid  surfaces,  the  wall  matching  procedure  of  Wilcox  was  employed  in  all  regions 
where  the  y+  of  the  first  cell  exceeds  ten.  At  high  speeds  the  flow  compressibility  effect  may  become  important  and 
it  can  lead  to  a  reduced  mixing  rate.  To  account  for  this  effect,  a  compressibility  correction  is  normally  introduced 
in  turbulence  models.  However,  all  present  compressibility  correction  models  seem  to  do  more  harm  than  good  in 
predictions  for  scramjet  engine  applications.  Therefore,  the  compressibility  correction  was  disabled  in  the  present 
simulation.  Also,  in  the  RANS  calculation,  the  turbulent  transport  of  momentum  is  modeled  by  the  eddy  viscosity 
model  while  the  turbulent  transport  of  energy  and  mass  is  modeled  using  the  gradient  diffusion  hypothesis.  The 
diffusion  rates  for  energy  and  mass  are  controlled  by  specifying  the  turbulent  Prandtl  (Prt)  and  Schmidt  (Sct) 
numbers.  Constant  values  for  these  numbers  are  usually  assumed  in  applications,  although  their  values  have  been 
found  to  vary  spatially.9’ 10  Calculations  performed  by  various  authors11’  12’  13  have  shown  strong  sensitivity  of  the 
solution  to  these  numbers,  especially  to  Sct.  Therefore,  extreme  care  should  be  taken  when  attempting  to 
characterize  scramjet  engines  with  constant  turbulent  transport  coefficients.  This  is  especially  true  for  the  present 
study  because  the  turbulent  mixing  of  fuel  and  air  is  one  of  the  most  important  phenomena  for  conditions  near 
blowout,  and  the  value  of  Sct  controls  the  extent  of  turbulent  mixing.  In  this  study,  Prt  and  Sct  were  chosen  to  be 
0.89  and  0.5,  respectively,  unless  otherwise  specified. 
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IV.  Results  and  Discussion 


A.  Validation  study 

To  demonstrate  the  accuracy  of  the  reduced  kinetic  mechanism  for  ethylene  and  also  the  correctness  of  its 
implementation  into  the  VULCAN  code,  a  validation  study  was  first  carried  out  by  employing  the  modified 
VULCAN  code  to  predict  the  ignition  delay  time  (Tign)  for  an  auto-ignition  problem  using  a  range  of  initial 


(a) 


(b) 


(c) 

Figure  2.  Comparisons  of  ignition  delay  time  at  the  equivalence  ratio  of  (a)  9=0.6,  (b)  9=1.0,  and  (c)  9=1.5. 
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temperatures,  pressures,  and  equivalence  ratios.  For  this  study,  only  a  single  grid  cell  was  needed  and  the  problem 
thus  had  constant  volume.  All  simulations  were  conducted  in  the  DOD  ASC  MSRC's  SGI  Origin  3900  machines. 
The  time  steps  in  the  transient  simulation  strategy  were  different  for  different  cases  and  they  varied  from  10  9to  10'6 
second.  Temperature  became  steady  between  5,000  and  50,000  steps  and  the  runtime  was  less  than  15  minutes. 
Figure  2  shows  the  comparisons  of  xign  between  the  15-step  reduced  and  detailed  kinetic  mechanism  calculations  for 
equivalence  ratios  of  (p=0.6,  1.0.,  and  1.5.  For  a  specific  cp,  a  total  of  nine  cases  were  considered  for  the  reduced 
kinetic  mechanism  with  initial  temperatures  of  1200  K,  1500  K,  and  1800  K  and  initial  pressures  of  0.1  atm,  1.0 
atm,  and  30.0  atm.  It  is  seen  that  the  Tign  from  the  reduced  kinetic  mechanism  agrees  very  closely  with  that  from  the 
detailed  one  for  the  entire  parameter  range. 

B.  AFRL/PRA  cavity  combustor  simulation  with  the  15-step  reduced  kinetic  mechanism 

AFRL/PRA  direct-connect  cavity  combustor  test  rig  features  a  facility  nozzle,  an  isolator  section,  a  cavity 
section,  and  a  post-cavity  section.  Figure  3  shows  half  of  the  geometry  configuration  of  this  test  rig.  A  close-up 
view  of  the  cavity  section  is  provided  in  Fig.  4.  Room- temperature  ethylene  in  the  test  was  fed  into  the  cavity  from 
ten  injectors  placed  on  the  cavity  ramp  wall.  The  air  mass  flowrate  was  varied  by  adjusting  the  plenum  pressure 
upstream  of  the  facility  nozzle.  At  each  test  condition,  the  fuel  flowrate  was  gradually  increased  until  the  cavity 
was  ignited.  Once  the  cavity  was  lit,  the  fuel  flow  was  gradually  decreased  to  the  point  where  the  cavity  blew  out. 
For  the  current  numerical  study,  the  selected  inflow  test  condition  was  a  Mach  2  flow  with  a  total  pressure  and 
temperature  of  414  kPa  and  589  K,  respectively.  The  incoming  air  mass  flowrate  was  3.1  kg/sec  and  the  measured 
fuel  flowrate  at  lean  blowout  was  0.46  kg/sec.  This  test  datum  will  be  used  to  validate  the  present  numerical  results. 


Post  Cavity 


Figure  3.  Geometric  configuration  of  AFRL/PRA  cavity  combustor. 


Figure  4.  Close-up  view  of  the  cavity  section. 
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In  an  effort  to  reduce  the  computational  time  required  to  perform  the  simulations  without  sacrificing  grid 
resolution,  the  existence  of  symmetry  planes  between  adjacent  fuel  injectors  was  presumed.  This  approximation  is 
rigorously  correct  in  the  limit  of  an  infinite  array  of  injectors  where  the  combustor  side  walls  are  non-existent.  With 
this  presumption,  only  a  slice  of  the  combustor  facility  with  a  half  injector  was  needed  in  the  modeling.  Figure  5 
shows  the  computational  domain  used  for  the  AFRL/PRA  combustor  simulation.  The  total  grid  size  is  about  one 
million.  Figure  6  shows  a  close-up  view  of  the  grid  in  the  cavity  region  with  the  injector  tube.  In  order  to  check  the 
performance  of  the  reduced  kinetic  model  and  also  examine  the  flame  stability  limits,  simulations  with  different  fuel 
flowrates  were  carried  out  in  this  study.  The  lean  blowout  case  with  a  fuel  flowrate  of  0.46  kg/sec  was  considered 
as  the  baseline  case.  Three  more  cases  with  fuel  flowrates  of  4,  5,  and  10  times  the  baseline  flowrate  were 
investigated  and  they  are  labeled  as  the  4X,  5X,  and  10X  cases,  respectively,  in  the  following  discussion.  These 
four  cases  with  different  fuel  flowrates  cover  combustion  from  the  fuel-lean  to  fuel-rich  conditions. 


Figure  5.  Computational  domain  for  the  ARFL/PRA  cavity  combustor. 


Figure  6.  Close-up  view  of  the  grid  in  the  cavity  section. 


AFRL/PRA  cavity  combustor  was  designed  to  study  lean  and  rich  flame  blowout  only  in  the  cavity  and  it  was 
not  intended  to  achieve  sustained  combustion  in  the  core  flow.  Therefore,  a  change  of  fuel  flowrate  has  an  impact 
only  on  the  flow  in  the  cavity  region.  Thus,  to  further  reduce  the  computational  cost,  the  computational  domain  in 
Fig.  5  was  broken  into  two  separate  regions,  with  one  containing  the  facility  nozzle  and  another  containing  the 
isolator  and  cavity  combustor.  The  simulation  in  the  facility  nozzle  domain  was  conducted  first.  The  flow  profile 
at  the  exit  plane  from  this  simulation  was  then  saved  and  used  as  the  inflow  condition  for  the  simulation  containing 
the  domain  comprised  of  isolator  and  cavity  combustor.  The  simulation  for  the  facility  nozzle  needs  to  be  redone 
only  when  different  air  inflow  conditions  are  used.  Due  to  the  still  relatively  large  computational  load,  16  parallel 
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processors  were  used  for  each  simulation  and  they  were  from  the  DOD  ASC  MSRC's  SGI  Origin  3900  or  Compaq 
Sc-45  machines.  For  a  run  with  combustion,  ignition  was  achieved  by  introducing  an  ignition  source  in  the  cavity 
region.  This  source  was  removed  after  several  hundred  flow  solver  iterations.  All  solutions  reached  steady  state 
and  their  mass  conservation  errors  were  controlled  below  0.1%.  Except  for  the  initial  stage  of  the  solution  process, 
the  CFL  number  was  usually  set  to  be  4  or  5.  For  the  cases  with  very  small  fuel  flowrates  in  the  injector,  the  CFF 
number  had  to  be  reduced  to  1  or  2. 

Figures  7,  8,  9,  and  10  show  the  local  equivalence  ratio,  temperature,  and  C02  mass  fraction  contours  at  the 
baseline,  4X,  5X,  and  10  X  cases,  respectively.  The  surfaces  where  the  contours  are  shown  include  the  cavity  wall, 
injector  wall,  and  the  symmetry  plane  across  the  injector.  To  help  in  the  result  analysis,  the  stoichiometric  lines  are 
also  plotted  for  each  contour  and  they  are  represented  as  the  white  lines.  For  the  baseline  case  (Fig.  7),  fuel  in  the 
injector  is  not  choked  and  it  is  injected  into  the  cavity  at  very  low  speed  (M-0.1).  Except  in  the  region  close  to  the 
injector  port,  the  cavity  is  filled  with  air  as  seen  in  Fig.  7(a).  Under  such  a  fuel-lean  condition,  the  present  study 
predicted  blowout,  which  is  consistent  with  the  experimental  observation.  The  temperature  in  the  cavity  is  found  to 
be  approximately  equal  to  the  total  temperature  (=  589  K)  of  the  incoming  air  flow.  Due  to  the  flameout,  the 
combustion  product  C02  is  not  present  in  the  cavity  region  (Fig.  7(c)). 
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Figure  7.  Contours  of  (a)  local  equivalence  ratio,  (b)  temperature,  and  (c)  C02  fraction  for  the  baseline  case. 

For  the  4X  case  (Fig.  8),  the  fuel  flowrate  is  increased  to  be  four  times  higher  and  the  fuel  is  recirculated  into 
most  of  the  cavity  region,  as  seen  in  Fig.  8(a).  The  cavity  flow  is  usually  characterized  by  a  shear  layer  and  a 
recirculation  zone.  On  the  symmetry  plane,  the  stoichiometric  line  lies  inside  the  shear  layer,  which  separates  the 
recirculation  zone  and  the  outer  freestream  of  air.  The  flow  near  the  cavity  wall  is  found  to  be  fuel  rich  (red  color). 
Due  to  the  extensive  mixing  of  fuel  and  air,  sustained  combustion  was  predicted  and  the  flame  was  stabilized  in  the 
cavity  center  region,  including  parts  of  the  shear  layer  and  the  recirculation  zone,  as  shown  in  Fig.  8(b).  In  addition 
to  temperature,  C02  may  also  be  used  as  an  indicator  of  combustion  intensity.  However,  under  fuel-rich  or  high 
flame  temperature  condition,  C02  dissociates  into  CO  and  its  amount  is  reduced.  As  a  result,  a  high  C02  mass 
fraction  is  observed  near  the  stoichiometric  line  instead  of  in  the  region  where  the  peak  flame  temperature  is  located 
(Fig.  8(c)). 

As  the  fuel  amount  is  further  increased  to  five  times  higher  than  that  of  the  baseline  case  (Fig.  9),  the  flow  in  the 
recirculation  zone  becomes  more  fuel  rich  and  it  pushes  the  stoichiometric  line  upward,  especially  in  the  ramped  aft 
wall  region,  as  seen  in  Fig.  9(a).  Consequently,  the  flame  region  also  moves  to  the  aft  wall  (Fig.  9(b)).  Due  to  the 
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intensified  fuel-rich  condition,  combustion  in  the  center  of  the  recirculation  zone  becomes  weak  and  the 
corresponding  flow  temperature  is  reduced.  As  in  the  4X  case,  CO2  is  produced  mainly  near  the  vicinity  of  the 
stoichiometric  line  with  a  local  equivalence  ratio  less  than  2  as  shown  in  Fig.  9(c). 
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Figure  9.  Contours  of  (a)  local  equivalence  ratio,  (b)  temperature,  and  (c)  C02  mass  fraction  for  the  5X  case. 
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Figure  10  provides  the  contours  for  the  fuel  flowrate  of  ten  times  higher  than  that  of  the  baseline  case.  At  this 
condition  most  of  the  cavity  is  flooded  with  fuel  and  the  stoichiometric  line  is  pushed  into  the  higher  velocity  region 
near  the  freestream  at  the  top  of  the  cavity  (Fig.  10(a)).  The  flame  remains  lit  but  it  is  contained  in  a  small  region 
near  the  ramped  aft  wall  and  surrounded  by  the  fuel-rich  flow  on  the  left  and  high-speed  freestream  on  the  top  as 
seen  in  Fig.  10(b).  It  is  expected  that  a  further  increase  of  the  fuel  flowrate  would  result  in  an  even  smaller  flame 
region  and  colder  recirculation  zone,  eventually  leading  to  fuel-rich  blowout.  The  clearly  diminishing  combustion 
for  the  10X  case  is  also  revealed  by  the  C02  contour  in  Fig.  10(c);  the  pronounced  C02  is  now  only  found  in  a 
region  near  the  upper  portion  of  the  aft  wall. 
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Figure  10.  Contours  of  (a)  local  equivalence  ratio,  (b)  temperature,  and  (c)  C02  mass  fraction  for  the  10X 

case. 

In  the  above  simulations,  Sct  was  set  at  0.5.  As  indicated  in  the  previous  section,  Sct  controls  the  turbulent 
mixing  of  fuel  and  air  and  its  value  is  expected  to  have  a  strong  impact  on  the  turbulent  combustion  as  well  as  the 
flame  blowout  prediction.  In  this  study,  this  impact  was  demonstrated  by  considering  three  different  Scts  for  the 
fuel  flowrate  corresponding  to  the  4X  case.  Figures  11,  12,  and  13  provide  comparisons  of  local  equivalence  ratio, 
temperature,  and  C02  mass  fraction  contours  for  Sctof  0.3,  0.5,  and  1.0.  The  extent  of  mixing  of  fuel  and  air  is  best 
expressed  by  the  local  equivalence  ratio  contour  (Fig.  11).  At  Sct=0.3  (Fig.  11(a)),  fuel  and  air  mixing  is  very 
extensive  and  the  stoichiometric  line  spans  the  center  region  of  the  cavity.  The  flow  becomes  fuel-rich  only  in  a 
region  very  close  to  the  cavity  floor.  As  Sct  is  increased  from  0.3  to  0.5,  it  seems  that  the  fuel  and  air  mixing  slows 
down  (Fig.  11(b)).  The  fuel-rich  region  near  the  cavity  floor  becomes  thicker  and  it  expands  to  the  front  wall  of  the 
cavity.  As  Sct  is  further  increased  to  1.0  (Fig.  11(c)),  the  mixing  of  fuel  and  air  becomes  less  efficient.  The  fuel-rich 
region  expands  to  the  cavity  center  and  pushes  the  stoichiometric  line  upward  to  near  the  freestream  at  the  top  of  the 
cavity.  From  Fig.  11,  the  mixing  extent  of  fuel  and  air  is  found  to  be  related  to  the  size  of  the  red  fuel-rich  region. 
The  larger  the  red  region,  the  poorer  the  mixing  of  fuel  and  air. 

Figure  12  shows  the  temperature  contours  for  three  different  Scts.  At  Sct=0.3,  the  enhanced  mixing  of  fuel  and 
air  results  in  strong  combustion  across  most  of  the  cavity  region,  as  seen  in  Fig.  12(a).  The  corresponding  flame 
temperature  reaches  above  2600  K.  As  Sct  is  increased  from  0.3  to  0.5,  combustion  becomes  weaker  and  flame 
region  size  shrinks,  due  to  the  reduced  mixing  of  fuel  and  air.  At  Sct=1.0  (12(c)),  the  front  half  of  the  cavity  is 
flooded  with  fuel  and  the  flow  temperature  in  the  recirculation  zone  is  reduced  significantly.  Although  the  flame  is 
sustained,  it  is  pushed  to  a  small  region  near  the  injector  and  the  peak  temperature  is  reduced  to  about  2000  K. 
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Figure  12.  Temperature  contours  for  Sct  numbers  of  (a)  0.3,  (b)  0.5,  and  (c)  1.0 
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Figure  13  shows  the  C02  mass  fraction  contours  for  three  different  Scts.  As  indicated  before,  a  high  flame 
temperature  implies  strong  dissociation  of  C02  into  CO.  As  a  result,  a  relatively  small  quantity  of  C02  is  found  in 
the  cavity  region  at  Sct=0.3  (Fig.  13(a)),  while  a  significant  quantity  of  C02  is  observed  at  Sct=1.0  in  Fig.  13(c). 
The  results  shown  in  the  these  figures  clearly  demonstrate  a  strong  sensitivity  of  combustion  to  Sct.  In  theory,  the 
Sct  for  a  specific  application  should  be  determined  based  on  the  calibration  of  experimental  data.  In  reality, 
however,  the  rule  of  thumb  is  to  set  Sct=0.5.  This  value  seems  to  work  well  for  many  applications  involving 
scramjet  engines. 
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Figure  13.  C02  mass  fraction  contours  for  Sct  numbers  of  (a)  0.3,  (b)  0.5,  and  (c)  1.0 


C.  Comparisons  of  results  among  the  reduced  mechanism  and  two  global  mechanisms 

Currently,  two  most  widely  used  ethylene  combustion  models  to  support  the  HyTech  program  are  the  10-step 
quasi-global  reaction  mechanism  developed  by  Singh  and  Jachimowski14  and  the  3-step  global  mechanism 
developed  by  researchers  at  AFRL.13  These  two  models  have  been  validated  over  a  certain  range  of  initial  pressure, 
temperature,  and  equivalence  ratio  conditions  based  on  the  Xign.  While  these  two  models  are  very  simple  to  apply, 
studies  at  AFRL  and  other  organizations  have  shown  great  discrepancies  between  experiment  measurements  and  the 
results  obtained  using  these  two  models.  To  assess  their  validity  in  the  present  work  on  flame  stability,  these 
models  were  applied  to  simulate  the  same  ARFL/PRA  cavity  combustor  as  that  for  the  15-step  reduced  mechanism. 
With  the  10-step  model,  sustained  combustion  was  not  achieved,  no  matter  what  fuel  flowrate  was  applied  and  what 
ignition  approach  was  employed.  With  the  3-step  model,  sustained  combustion  was  predicted  only  for  the  case  with 
the  fuel  flowrate  not  smaller  than  that  corresponding  to  the  4  X  case  and  the  Sct  number  not  greater  than  0.5.  This 
implies  that  the  flames  predicted  by  the  two  global  models  blow  out  more  easily  than  that  predicted  by  the  15-step 
reduced  model. 

To  understand  the  reason  behind  this  phenomenon,  the  extinction  residence  time  (xext)  in  a  PSR  was  calculated 
using  the  CHEMKIN  software  package  for  three  different  reaction  models.  The  initial  pressure  and  temperature  in 
the  PSR  were  chosen  to  be  close  to  those  in  the  core  flow  near  the  cavity  region,  and  they  were  0.5  atm  and  325  K, 
respectively.  Table  1  shows  the  comparison  of  the  xext  at  different  cps  for  three  reaction  models.  Overall  the  xext 
predicted  by  the  two  global  models  are  one  or  two  orders  of  magnitude  higher  than  that  of  the  15-step  reduced 
model,  which  implies  that  the  flame  calculated  from  the  15-step  reduced  model  is  more  difficult  to  blow  out,  due  to 
the  shorter  extinction  time.  The  large  discrepancy  in  the  Xext  from  different  models  should  perhaps  not  be  surprising 
because  the  two  global  models  were  developed  to  match  the  xign  only  within  a  certain  range  of  the  initial  conditions, 
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and  they  were  NEVER  intended  to  match  the  Xext  in  any  condition.  The  extinction  and  the  ignition  processes  are 
intrinsically  different,  in  that  a  radical  pool  already  exists  for  extinction,  but  it  has  to  be  built  up  in  ignition.  For 
predictions  of  blow-out,  the  flame  already  exists  and  extinction  time  is  of  more  relevance  than  ignition  time. 
Therefore,  the  two  global  models  are  inadequate  for  flame  stability  predictions. 

Table  1.  Comparison  of  xext  (sec.)  for  three  kinetic  mechanisms 


(p 

Kinetic  mechanism 

3-step 

10-step 

15-step 

0.5 

5.74E-02 

3.05E-02 

6.12E-04 

1.0 

2.87E-03 

3.37E-03 

7.19E-05 

2.0 

1.35E-03 

1.94E-03 

3.19E-04 

Regarding  the  reason  why  the  10- step  model  failed  to  predict  sustained  combustion,  the  speculation  is  that  the 
Text  from  the  10- step  mechanism  is  longer  than  those  of  other  models,  as  seen  in  Table  1  for  the  region  most 
favorable  for  flame  formation.  The  local  equivalence  ratio  in  this  region  is  normally  around  1.0  or  slightly  higher, 
as  shown  in  the  previous  figures.  Another  possible  reason  is  that  the  10-step  model  contains  several  additional 
radical  species  (OH,  H,  and  O)  than  the  3-step  model.  As  a  result,  the  heat  release  and  temperature  predicted  by  the 
10- step  model  are  usually  lower  and  flameholding  is  more  difficult. 

Even  when  sustained  combustion  is  achieved,  the  flame  structures  predicted  by  different  models  may  be 
different.  Figures  14  and  15  show  comparisons  of  local  equivalence  ratio  and  temperature  contours  between  the 
reduced  and  3-step  models  for  Sct=0.5  and  the  fuel  flowrate  corresponding  to  the  5X  case.  While  the  local 
equivalence  ratio  contours  (Fig.  14)  predicted  by  the  two  models  are  very  similar,  their  temperature  contours  (Fig. 
15)  are  quite  different.  The  flame  predicted  by  the  3-step  model  is  seen  to  extend  much  closer  to  the  cavity  front 
wall  than  that  from  the  reduced  model.  In  addition,  the  flame  peak  temperature  predicted  by  the  3 -step  model  is 
also  higher.  One  major  reason  for  this  is  the  absence  of  radical  species  in  the  3-step  model.  The  energy  consumed 
in  radical  formation  is  higher  than  the  energy  released  in  radical  combination.  Neglecting  radical  results  in  an  over¬ 
prediction  of  the  overall  energy  release  in  the  flow  and  flame  temperature.  In  addition,  the  3-step  model  is  only 
valid  over  a  narrow  range  of  conditions.  Calculations  made  outside  of  these  validation  conditions  may  incur 
significant  errors  in  the  results.  For  example,  CHEMKIN  was  used  to  calculate  the  Xignover  a  range  of  pressures  and 
equivalence  ratios.  Both  the  3-step  and  10-step  models  predicted  Xignthat  are  at  least  an  order  of  magnitude  shorter 
than  that  from  the  detailed  kinetic  mechanism,  once  the  temperature  dropped  below  1200  K.  A  short  xign  implies 
enhanced  combustion.  Based  on  the  above  analysis,  the  15-step  reduced  kinetic  mechanism  is  preferred  for 
ethylene  combustion  simulations  because  this  model  has  been  validated  for  many  parameters  such  as  the  xign,  xext, 
laminar  flame  speed,  etc.  over  a  wide  range  of  flow  conditions. 
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Figure  15.  Temperature  contours  for  (a)  the  reduced  kinetic  mechanism  and  (b)  3-step  global  mechanism. 


V.  Conclusions 

The  VULCAN  code  along  with  a  reduced  chemical  kinetic  mechanism  for  ethylene  oxidation  has  been  applied 
to  model  the  AFRL/PRA  cavity  combustor  and  predict  cavity- stabilized  flames  at  selected  test  run  conditions. 
Several  fuel  flowrates  have  been  considered,  and  combustion  in  the  cavity  involves  fuel-lean  as  well  as  fuel-rich 
conditions.  The  predicted  fuel  flowrate  at  the  lean  blowout  condition  seems  to  match  the  experimental  data.  By 
considering  several  cases  with  different  Scts,  the  cavity  combustion  has  been  found  to  be  sensitive  to  Sct.  A  value  of 
0.5  for  Sct  appears  to  be  adequate  for  the  present  cavity  combustor.  To  demonstrate  the  importance  of  using  the 
reduced  kinetic  mechanism,  two  widely-used  global  kinetic  models  have  been  applied  to  simulate  the  present  cavity 
combustor  and  they  have  been  found  to  be  inadequate  for  applications  involving  flame  stability  prediction.  Because 
the  15-step  reduced  kinetic  mechanism  has  been  validated  for  many  parameters  over  a  wide  range  of  flow 
conditions,  it  is  preferred  for  ethylene  combustion  simulations. 
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